Scaling law for the transient behavior of type-II neuron models 
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We study the transient regime of type-II biophysical neuron models and determine the scaling 
behavior of relaxation times r near but below the repetitive firing critical current, r ~ C{Ic — 
7)"'^. For both the Hodgkin- Huxley and Morris-Lecar models we find that the critical exponent 
is independent of the numerical integration time step and that both systems belong to the same 
universality class, with A — 1/2. For appropriately chosen parameters, the FitzHugh-Nagumo 
model presents the same generic transient behavior, but the critical region is significantly smaller. 
We propose an experiment that may reveal nontrivial critical exponents in the squid axon. 
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I. INTRODUCTION 

The search for biophysical models for information pro- 
cessing systems has led to a variety of model neurons 
which describe the dynamics of the membrane poten- 
tial. Collective properties arise from their interaction 
through several possible architectures and types of cou- 
plings. These can be chemical, voltage-gated synapses, 
simpler proteic electrical connections, or even just elec- 
trical cphaptic interactions arising between neighboring 
nerve fibers. Maybe the most striking feature of a neu- 
ron is the threshold of the stimulus that separates spik- 
ing from nonspiking regimes. Spiking neurons generate 
taletale signatures which have served as the basis for 
frequency-dependent neural codes, an idea that can be 
traced back to the work of Adrian in the 1920s Al- 
though of paramount importance to neural dynamics, 
spike frequencies do not tell the complete story. Sub- 
tle computations may arise from subthreshold dynamics 
such as for examples in the early stages of the mammalian 
visual system, olfactory bulb, and cortex. In many cases 
the key to the information dynamics lies in the transients, 
either below the current threshold to generate action po- 
tentials or the threshold to generate infinite sequences of 
spikes. In this paper we investigate transient spike trains 
of single model neurons since this might have a bearing 
on the collective behavior, i.e., computational capabili- 
ties, of subthreshold assemblies of neurons. 

A dynamical system approach reveals that universal 
bifurcation scenarios for the firing behavior appear irre- 
spective of the specific membrane ion channels and mi- 
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croscopic details involved. Full characterization of these 
bifurcation routes is important for a deeper understand- 
ing of how they affect the firing behavior and possible 
implementation of neural codes. For example, it is now 
acknowledged that bistable behavior and the small range 
of firing frequencies in neurons that undergo subcritical 
Hopf bifurcations prevent a robust use of a pure fre- 
quency code. 

The transient behavior of neuron models has not re- 
ceived much attention in either experimental or theoret- 
ical studies. In this paper we study the divergence of 
transient times in a class of conductance-based models 
and show that it follows a universal critical behavior. We 
propose an experiment that may test our theoretical pre- 
dictions and discuss how neurons could employ transients 
for computational purposes. 



II. TRANSIENTS IN TYPE-II MODELS 

The Hodgkin-Huxley (HH) model is a biophysically 
motivated system of four coupled nonlinear differential 
equations that describe the dynamics of the membrane 
potential V of the squid giant axon (e.g., 0, [1]): 

= GNam^h{ENa'V) + GKnHEK -V) 

at 

+Gl {EL-V) + I{t), 
^ = a,^{V){l-x,)-P,^{V)x, . (1) 

Xi stand for the three gate variables Xi = m, h, and n 
describing the activation of ionic channels and , Pxi 
are voltage-dependent transition rates Q : 
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where voltages are expressed in mV, rates in ms ^, 
C = 1 /^F/cm^ is the specific membrane capacitance, 
Gjva = 120 mS/cm2, Gk = 36 mS/cm^, and Gl = 
0.3 mS/cm^ are ionic conductances per unit area, and 
Ejsia = 115 mV, Ek = -12 mV, and El = 10.6 mV arc 
reversal potentials. 

The HH system plays a fundamental role in the field 
of neurophysiology and computational ncuroscicncc since 
it defines the class of conductance-based models. As a 
function of a constant applied current /, the HH model 
undergoes a subcritical Hopf bifurcation at I = Ih, above 
which the fixed point solution (membrane potential at 
rest) is no longer stable and trajectories are attracted to 
a stable limit cycle, leading to repetitive firing (infinite 
train of action potentials). 

The sudden jump to a periodic behavior with nonzero 
frequency / is analogous to a first-order phase transition 
and is referred to as type-II behavior in the neuroscience 
literature As in first-order phase transitions, coex- 
istence also appears in type-II behavior. Just below the 
Hopf bifurcation, the stable fixed point coexists with a 
stable limit cycle, their basins of attraction being sepa- 
rated by an unstable limit cycle Q . Both limit cycles are 
created in a saddle-node (or "fold" ) bifurcation of cycles 
at I = Ic < Ih ■ In our analogy with equilibrium phase 
transitions, Ic would correspond to a spinodal point. If 
the fixed point at / = is perturbed by the application 
of a constant current near but below Ic, several spikes 
may appear before the system returns to the new resting 
state (Fig. [I]). 

In the original version of the Morris-Lecar (ML) 
model , a system with two coupled nonlinear ordinary 
differential equaions (ODEs) used to describe action po- 
tentials in a barnacle motor fiber, the relevant bifurca- 
tion at the onset of repetitive firing is a saddle-node one. 
That means that the spiking frequency varies continu- 
ously from Ic as f (X (/ — Ic)^, with /? ~ 1/2, which 
is similar to a mean field second-order phase transition 
behavior if we think of / as the order parameter. This 
transition is called type-I behavior in the neuroscience 
literature and does not present a slow transient phe- 
nomenon similar to Fig. 1. However, the Morris-Lecar 
system has also been used to describe cells which present 
a type-II behavior 0, Q , which occurs for the equations 
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FIG. 1; Examples of transient behavior near Ic for the HH 
model. The constant step current is applied at t = 10 ms. 
The estimated critical current is Ic = 6.26422125685 fiA/crn^ 
for an integration time step dt — 0.01 ms. 
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when, for example, the values of the parameters are 
chosen as Gca = 1-1 mS/cm^, Gk = 2.0 mS/cm^, 
Gl = 0.5 mS/cm2, Eca = 100 mV, Ek = -70 mV, 
and El = — 50 mV. Then, large transient times are also 
observed (Fig. [2]). 

The FitzHugh-Nagumo (FHN) system 



dV 
'dt 

dw 



V{V -a){l-V)-w + I, 



(4) 



has been proposed as a low-dimensional toy model that 
represents the type-II behavior of the HH and other ex- 
citable systems. We verified that the transient behavior 
here reported is not seen with the usual parameters [1]. 
However, the FHN model can reproduce the HH tran- 
sient behavior if one chooses parameters near a = 0.5, 
7 = 4.2, and e = 0.01 (Fig.©. 

In this paper we show that the long relaxation times 
in type-II models are a consequence of the changes in 
phase space which occur near the creation of the limit 
cycles. Moreover, this leads to a scaling relation whose 
exponent can be predicted, in agreement with numerical 
simulations. 
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FIG. 2: Examples of transient behavior near Ic for the 
ML system. The constant step current is applied at 
t = IQ ms. The estimated critical current is Ic ~ 
24.84134676279 ^K/cvc? for an integration time step dt — 
0.01 ms. 
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FIG. 4: Relaxation times for the HH model as a function of 
the distance to critical current for different integration times: 
Tmax — 10'' ms (open circles) and Tmax — 10^ ms (filled 
circles) . 
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FIG. 3: Examples of transient behavior near Ic for the FHN 
system. The constant step current is applied a.t t = 10. The 
estimated critical current is Ic = 0.1025447183127 for an inte- 
gration time step dt — 0.01 (all quantities in arbitrary units). 
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FIG. 5: Relaxation times for the Morris-Lecar model as a 
function of the distance to critical current for different inte- 
gration times: dt = 0.01 ms (filled circles) and dt — 0.1 ms 
(open circles). 



III. SCALING LAW 

The relaxation time r may be defined as the time until 
the last spike, or the time until the membrane voltage 
stays within a small distance from the resting potential 
(these two times arc very similar near Ic). When we plot 
T as a function of Ic-- 1 we find a power law divergence of 
the relaxation time, t ~ C{Ic — I)^'^ , where A is similar 
to a dynamic critical exponent. The A exponent charac- 
terizes the "critical slowing down" behavior near the bi- 
furcation. We expect that A is a universal exponent but 
we are not aware that this exponent has been measured 
for neuron models. Here we report the measured expo- 



nents for these three type-II biophysical neuron models, 
finding very good agreement between them. 

We integrated the equations using a standard fourth- 
order Runge-Kutta algorithm and determined r by mea- 
suring the time interval from the onset of the current step 
to a near stop of the flow [|x| < 10~^, where x is the ve- 
locity vector in phase space: x = (w, v) for the ML and 
FHN systems, and x = {wjjh, h,n) for the HH model]. 
As opposed to the Hopf bifurcation, the fold bifurcation 
cannot be obtained analytically, so Ic was estimated nu- 
merically after integration of the ODEs up to a (long) 
maximum time Tmax- The determination of the criti- 
cal current is sensitive to Tmax, but in practice this only 
limits the range of validity of the power law (see Fig.d]). 
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FIG. 6: Relaxation times for the FitzHugh-Nagumo model 
as a function of the distance to critical current. Note that 
the power law becomes visible only very close to the fold 
bifurcation {h — I ^ 10~^). 



We have employed Tmax ~ 10^ ms and dt ~ 0.01 ms, 
unless otherwise stated. The estimated eritieal currents 
Ic{Tmax,dt) quoted in the figure captions are very pre- 
cisely determined given T„iax and the integration step 
size dt. 

The critical exponent is determined from the plot r vs 
Ic-I- We found A = 0.47 for the HH system (Fig.gl) and 
A — 0.49 for the ML model (Fig. O, irrespective of the 
size of the integration step dt. This suggests a universal 
exponent A = 1/2. Obtaining long transients in the 
FHN model has proved numerically more difficult, since 
the phenomenon occurs only very close to Ic (Fig. (H]). 
Nonetheless, we have obtained the exponent A = 0.48. 

Figure [Tja) shows that for /c an unstable 

limit cycle and a stable one coexist and surround a sta- 
ble fixed point. The fixed point for / = lies outside 
both limit cycles [Fig. [7jb)]. Therefore, when the cur- 
rent is abruptly changed to / < /c, the fixed point is 
displaced to a region within a ghost limit cycle, and the 
transient is completely dominated by the time it takes 
for the system to overcome it. The ghost is a natural 
consequence of the system being immediately below a 
saddle-node bifurcation of cycles, and can be character- 
ized by the vanishingly small flow component normal to 
the half-stable limit cycle that is created at I = Ic- It ef- 
fectively works as a one-dimensional extended bottleneck 
through which the system must pass before reaching the 
fixed point (see Fig. [8] for a caricature). If one considers 
an analogous system in polar coordinates 0] 9 = f{r,9), 
r = iir + r'^, where /(r, 9) > ,V7' > 0, it is clear 
that for /i < /ic = — 1/4 the time for the system to over- 
come the ghost at r = 1/2 scales as t ~ {fic — m)^^^^ 
(the Hopf bifurcation occurring only at = 0). Solu- 
tions for the transient behavior have been obtained by 
Tonnelier Q for McKean's piecewise linear version of the 
FHN model However, the scaling behavior has not 
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FIG. 7: Phase portraits of the ML model for I slightly above 
(a) and below (b) Ic- The long transient is dominated by 
the time it takes to pass through the region where the limit 
cycles are about to emerge. SPO (UPO) = Stable (Unstable) 
Periodic Orbit, SFP = Stable Fixed Point. Horizontal axes 
(V) in mV and vertical axes (W) are dimensionless. 



been observed because the model has been studied in the 
absence of an external current. 

It should be clear that the current step is just a sim- 
ple way of putting the system outside the ghost limit 
cycle, but the scaling law for the transient is not re- 
stricted to this somewhat artificial protocol (even though 
it is very common in both theory and experiments). 
Close to the fold transition, any short-lived perturba- 
tion that is strong enough to make the system cross 
the ghost limit cycle will give rise to long transients 
back to the fixed point. We exemplify this with a bi- 
ologically plausible example in the HH model. Suppose 
the system is somehow maintained close to criticality at 
/ < /c (this could be achieved by several different pos- 
sible mechanisms, so we just fix /). In addition, assume 
the neuron undergoes a fast excitatory postsynaptic po- 
tential (EPSP) simulated by an injected synaptic current 

Isyn{t) = gsyn{t){ENa — V): 
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FIG. 8: Schematic bifurcation diagram for type II neuron 
models. The fold bifurcation occurs at 7c, while the Hopf 
bifurcation occurs at Ih- Owing to the onset of the current 
step (dot-dashed arrow), the fixed point for 7 = becomes 
the initial condition in a new phase portrait. The transient 
T (solid arrow) to reach the new fixed point is governed by 
the ghost limit cycle where the UPO and the SPO annihilate 
each other (see Fig. ^ . 



^{t'){9mt' /t^) cxp(— f'/rs), where is the Heaviside func- 
tion, g,n ~ 1 mS/cm^, = 0.5 ms, and t' = t — Iepsp^ 
where Iepsp = 50 ms is the time the EPSP is initiated. 
Results are shown in Fig. [51 whose similarity with Fig. [1] 
attests to the robustness of the effect (notice however 
that, differently from Fig. [1] in Fig. [5] the resting mem- 
brane potential is at the / 7^ fixed point before the 
perturbation). This opens interesting possibilities from 
the point of view of neuronal computation: the length of 
the transient response of a neuron "probed" by an EPSP 
could code for its internal state of excitability, that is, 
for how close it is to Ic- Naturally, this transient cod- 
ing would work only if the system is close to the fold 
bifurcation. We therefore could have another example in 
neuroscience of optimal information processing at criti- 
cality [HI, [H, [H, M [H, [H, ■ 

It is interesting to point out that this bottleneck effect 
is analogous to what occurs for typc-I neurons above the 
saddle-nodc bifurcation that leads to repetitive firing. In 
that case, however, the ghost results from the anihilation 
of fixed points (not limit cycles) and the period T of the 
limit cycles diverges as (/— /c)~^^^. This provides a com- 
plementary scenario connecting both classes of neurons: 
the transient of type-II models below Ic diverges with the 
same exponent as the period of type-I models above Ic, 
that is, A = /3. 



IV. CONCLUDING REMARKS 
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FIG. 9: Long transients appear if the system initially at rest 
with 7 < 7c is perturbed by an additional EPSP (see text 
for details). Solid lines are membrane potentials, while the 
dashed line is the synaptic conductance. 



dt 



GNam^h {ENa - V) + Gku" [Ek - V) 

+GL{EL-V)+I + hyn{t) . (5) 



The fast change in the synaptic conductance is 
given by Rail's alpha function [l^]: g{t') = 



We were unable to find a description of this scaling law 
behavior for transients in neuron models or the associated 
dynamic critical exponent in the literature. Some kind 
of critical slowing down for subthreshold oscillations has 
been reported experimentally in the squid axon [Tsj , but 
these authors examined the vicinity of a parametric sub- 
critical Hopf bifurcation, not a saddle-node bifurcation of 
cycles induced by external currents. We propose that a 
similar experiment with high-precision injected currents 
near could be used to check the power law found in 
the computational model. Since the area of the giant 
squid axon is of order ~ 1 cni^, to examine the critical 
regime requires that current fluctuations should be less 
than ~ 10'^ pA (see Fig. [3]). Even if the full critical regime 
seems to be hard to achieve, the initial divergence in the 
transient lifetime may provide an experimental check of 
our predictions. 

We emphasize that both in standard experiments and 
in our single-compartment model space clamping is used. 
It might happen though, as in spin systems, that for 
the extended real system without voltage clamp, or for a 
compartmental model with a large number of compart- 
ments, the exponents may differ from the values here 
reported, changing the universality class to one not de- 
scribed by a "mean field" approach. Interestingly, this 
would mean that collective properties within a nontriv- 
ial universality class could be observed at the level of a 
single axon. So whether this result can be verified exper- 
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imcntally hinges on the effects that spatially extended 
neurons may have on the robustness of this picture. 

Furthermore, noise could always play a role. In stud- 
ies of type-I intermittency in simple maps, the length of 
the "laminar phase" (l) in a chaotic regime is analogous 
to the transient in this work and diverges as (I) ^ iT^I'^ 
because of a zero-dimensional bottleneck as the distance 
e to a tangent bifurcation tends to zero 1£ . [20II . In the 
presence of additive noise with amplitude g 2l| , the scal- 
ing changes to (Q ~ where / is a universal 
function (see also [24| for recent extensions). It is con- 
ceivable that similar scaling relations could be obtained 
for T provided that the chaotic phase of intermittency 
theory could be replaced by some mechanism of "rein- 



jection" in typc-II neuron models. For instance, in the 
phenomenon of coherence resonance [23l | noise itself plays 
this reinjecting role. However, the excitable systems em- 
ployed are usually not close enough to the fold transition 
to exhibit long transients, so it would be interesting to 
investigate the effects of the scaling laws we report here 
in the resonance curves. These theoretical issues should 
be dealt with before engaging in an experimental search. 
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